/*  sdfsqrt.cpp  by K.Tsuru  */
// function ID 3007  DRADIX
/***********************************
SDouble class
It provides the square root sqrt(x).
************************************/
#ifndef SN_H
#include "sn.h"
#endif
static const char* func = "Sqrt";
SDouble Sqrt(const SDouble& x){
    if( !x.Sign(3007) ) return 0.0;
    if( x.Sign() < 0 ) x.SetError(x.DOMAIN_ERR, func, 3007);

    SDouble y, rs, d;
//It checks whether an exact value can be obtained by double or not.
//It costs a little time bacause the "dsqrt" is a short number.
    uint xfig = x.Last() - x.First();
    if(xfig < 4){
        SDouble dsqrt;
        y.SNError();    //reset error
        dsqrt = sqrt( doubleD(x, 0) );
        if(y.SNError()==y.NO_ERR){  //detect an error in doubleD()
            y = x - dsqrt*dsqrt;
            if(y.Sign(3007)== 0) return dsqrt;
        }
    }
// sqrt(x) = x*(1/sqrt(x))
    rs = RecSqrt(x);    // = 1/sqrt(x)
    y = x*rs;
/***************************************************************
It adds a correction.In some cases it rases up the precision
by about ten percent.It will be confirmed by comparing the two values
Sqrt(Pi()) in the effective figures 120 and 130.
 y' = sqrt(x-d) = sqrt(x)*sqrt(1-d/x)
 sqrt(x) = y'/sqrt(1-d/x) ~= y'(1+0.5*d/x) = y' + 0.5*d*y'/x
 ~= y'+0.5*d*(1/sqrt(x)) = y' + 0.5*d*rs
*****************************************************************/
#if 1
    if(y.PreferSpeed() == OFF){
        d = x - y*y;    // "d" has a few figures.
        if(d.Sign()) y += DsDiv(d, 2)*rs;
    }
#endif
    if(y.Verify()){     // check |x - y^2| << x ?
        d = x - y*y;
        if(!d.IsMLT(x)){  //It sees the relative error. d<<x?
            y.SetError(y.VERIFY, func, 3007);
        }
    }
    return y;
}